library(ggplot2)
  library(reshape2)
  library(dplyr)
  library(tidyr)
  library(GGally)
  library(grid)
  "%&%" = function(a,b) paste(a,b,sep="")
  source('/Volumes/im-lab/nas40t2/hwheeler/PrediXcan_CV/GTEx_2014-06013_release/transfers/PrediXmod/Paper_plots/multiplot.R')
  my.dir <- '/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/'
  fig.dir <- '~/GitHub/GenArch/GenArchPaper/Figures/'
  rna.dir <- my.dir %&% "gtex-rnaseq/"
  out.dir <- rna.dir %&% "ind-tissues-RPKM/"
  my.vol <- '/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/BSLMM_exp/'
  otd.dir <- '/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/gtex-OTD-CV-R2/'

Fig S1

Tissue-Wide joint global estimates

tislist <- scan(my.dir %&% 'rmTW.ten.tissue.list',sep="\n",what="character")[2:10]##rm cross-tissue from plot
ts <- data.frame()
for(tis in tislist){
  data <- read.table(my.dir %&% 'gtex-h2-estimates/GTEx.TW.' %&% tis %&% '.h2.all.models_FHSfdr0.05.Chr1-22_globalOtherChr.2015-10-06.txt',header=T,sep="\t") 
  explist <- scan(out.dir %&% tis %&% ".meanRPKMgt0.1_3samplesRPKMgt0_genelist","c")
  data <- dplyr::filter(data,ensid %in% explist)
  glo.jt <- select(data,tissue,glo.jt.h2,glo.jt.se) %>% mutate(glo.jt.h2=ifelse(is.na(glo.jt.h2),0,glo.jt.h2), glo.jt.se=ifelse(is.na(glo.jt.se),sample(glo.jt.se[is.na(glo.jt.se)==FALSE],replace=TRUE,size=length(glo.jt.se[is.na(glo.jt.se)==TRUE])),glo.jt.se)) %>% mutate(ymin = pmax(0, glo.jt.h2 - 2 * glo.jt.se), ymax = pmin(1, glo.jt.h2 + 2 * glo.jt.se)) %>% arrange(glo.jt.h2) %>% mutate(`CI > 0`=ymin>0,place=1:length(data$tissue)) 
  ts <- rbind(ts,glo.jt)
}
p<-ggplot(ts,aes(x=place,y=glo.jt.h2,ymin=ymin,ymax=ymax,color=`CI > 0`) ) + facet_wrap(~tissue,ncol=3) + geom_pointrange(col='gray')+geom_point()+ylab(expression("global h"^2)) + xlab(expression("genes ordered by joint global h"^2))+theme_bw() + coord_cartesian(ylim=c(0,1))

###calc % nonzero for each tissue
### ADD to plot legend
pvec<-vector()
h2vec<-vector()
a<-ts %>% select(tissue,`CI > 0`) %>% spread(tissue,`CI > 0`)
for(i in 1:9){
  tis<-colnames(a)[i]
  cat("\n\n---",tis,"---\n")
  print(table(a[,i]))
  per <- signif(table(a[,i])/sum(table(a[,i])),3)*100
  print(per)
  pvec <- c(pvec,per[2])
}
## 
## 
## --- Adipose-Subcutaneous ---
## 
## FALSE  TRUE 
## 14193    12 
## 
##   FALSE    TRUE 
## 99.9000  0.0845 
## 
## 
## --- Artery-Tibial ---
## 
## FALSE  TRUE 
## 13501     3 
## 
##    FALSE     TRUE 
## 100.0000   0.0222 
## 
## 
## --- Heart-LeftVentricle ---
## 
## FALSE 
## 13321 
## 
## FALSE 
##   100 
## 
## 
## --- Lung ---
## 
## FALSE  TRUE 
## 14773     2 
## 
##    FALSE     TRUE 
## 100.0000   0.0135 
## 
## 
## --- Muscle-Skeletal ---
## 
## FALSE  TRUE 
## 12738    95 
## 
## FALSE  TRUE 
## 99.30  0.74 
## 
## 
## --- Nerve-Tibial ---
## 
## FALSE  TRUE 
## 14509     1 
## 
##    FALSE     TRUE 
## 1.00e+02 6.89e-03 
## 
## 
## --- Skin-SunExposed(Lowerleg) ---
## 
## FALSE  TRUE 
## 14604    21 
## 
##  FALSE   TRUE 
## 99.900  0.144 
## 
## 
## --- Thyroid ---
## 
## FALSE  TRUE 
## 14640     2 
## 
##    FALSE     TRUE 
## 100.0000   0.0137 
## 
## 
## --- WholeBlood ---
## 
## FALSE  TRUE 
## 12067    93 
## 
##  FALSE   TRUE 
## 99.200  0.765
###calc mean h2 for each tissue
a<-ts %>% select(tissue,glo.jt.h2) %>% spread(tissue,glo.jt.h2)
for(i in 1:9){
  tis<-colnames(a)[i]
  meanh2 <- signif(mean(a[,i],na.rm=TRUE),3)
  cat(tis,"mean h2:",meanh2,"\n")
  h2vec <- c(h2vec,meanh2)
}
## Adipose-Subcutaneous mean h2: 0.0487 
## Artery-Tibial mean h2: 0.0546 
## Heart-LeftVentricle mean h2: 0.0431 
## Lung mean h2: 0.0742 
## Muscle-Skeletal mean h2: 0.0425 
## Nerve-Tibial mean h2: 0.0472 
## Skin-SunExposed(Lowerleg) mean h2: 0.046 
## Thyroid mean h2: 0.046 
## WholeBlood mean h2: 0.0498
pvec<-ifelse(is.na(pvec),0,pvec)
ann_text <- data.frame( glo.jt.h2 = rep(0.78,9), place = rep(5000,9), percent= pvec, mean_h2 = h2vec, tissue = factor(colnames(a)), ymin=rep(0.9,9), ymax=rep(0.9,9), `CI > 0`=rep(NA,9), se=rep(0.9,9))
p2<-p+geom_text(data=ann_text,aes(label=paste("mean_h^2 ==",mean_h2,sep="")),color="black",show_guide=F,parse=T,hjust=0,size=3) 
ann_text <- data.frame( glo.jt.h2 = rep(0.9,9), place = rep(5000,9), percent= pvec, mean_h2 = h2vec, tissue = factor(colnames(a)), ymin=rep(0.9,9), ymax=rep(0.9,9), `CI > 0`=rep(NA,9), se=rep(0.9,9))
p3<-p2+geom_text(data=ann_text,aes(label=paste("percent_TRUE ==",percent,sep="")),color="black",show_guide=F,parse=T,hjust=0,size=3)+ theme(legend.justification=c(0,1), legend.position=c(0,1))

p3

png(filename=fig.dir %&% "Fig-GTEx_TW_glo.jt.h2.png",width=720,height=480)
p3
dev.off()
## quartz_off_screen 
##                 2
tiff(filename=fig.dir %&% "Fig-GTEx_TW_glo.jt.h2.tiff",width=720,height=480)
p3
dev.off()
## quartz_off_screen 
##                 2

Fig S2

Plot GTEx Tissue-Wide Elastic Net results

tislist <- scan('/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/nine.tissue.list','c')
finalgdata <- data.frame()
for(tis in tislist){
  alpha1 <- read.table(my.dir %&% 'gtex-OTD-CV-R2/TW_' %&% tis %&% '_exp_10-foldCV_elasticNet_alpha1_hapmapSnpsCEU_all_chr1-22_2015-09-10.txt',header=TRUE) %>%  mutate(`1`=R2) %>% select(gene,`1`)
  ngenesall <- length(unique(alpha1$gene))
  alpha95 <- read.table(my.dir %&% 'gtex-OTD-CV-R2/TW_' %&% tis %&% '_exp_10-foldCV_elasticNet_alpha0.95_hapmapSnpsCEU_all_chr1-22_2015-09-10.txt',header=TRUE) %>% mutate(`0.95`=R2) %>% select(gene,`0.95`)
  alpha50 <- read.table(my.dir %&% 'gtex-OTD-CV-R2/TW_' %&% tis %&% '_exp_10-foldCV_elasticNet_alpha0.5_hapmapSnpsCEU_all_chr1-22_2015-09-10.txt',header=TRUE) %>% mutate(`0.50`=R2) %>% select(gene,`0.50`)
  alpha05 <- read.table(my.dir %&% 'gtex-OTD-CV-R2/TW_' %&% tis %&% '_exp_10-foldCV_elasticNet_alpha0.05_hapmapSnpsCEU_all_chr1-22_2015-09-10.txt',header=TRUE) %>% mutate(`0.05`=R2) %>% select(gene,`0.05`)
  data <- inner_join(alpha05,alpha50,by='gene')
  data <- inner_join(data,alpha95,by='gene')
  data <- inner_join(data,alpha1,by='gene')
  gdata <- gather(data,alpha,R2,2:4) %>% mutate(tissue=tis)
  finalgdata <- rbind(finalgdata,gdata)
}

p<-ggplot(finalgdata, aes(y = `1` - R2, x = `1`, group=alpha, color=alpha)) + facet_wrap(~tissue) + geom_point(show_guide = TRUE) + ylab(expression(paste("R"^2, " difference (LASSO - alpha)"))) + xlab(expression(paste("R"^2, " (LASSO)"))) +theme_bw(15)+ theme(legend.justification=c(0,1), legend.position=c(0,1))
p
## Warning: Removed 24133 rows containing missing values (geom_point).
## Warning: Removed 24510 rows containing missing values (geom_point).
## Warning: Removed 27268 rows containing missing values (geom_point).
## Warning: Removed 25140 rows containing missing values (geom_point).
## Warning: Removed 25335 rows containing missing values (geom_point).
## Warning: Removed 22325 rows containing missing values (geom_point).
## Warning: Removed 23495 rows containing missing values (geom_point).
## Warning: Removed 23178 rows containing missing values (geom_point).
## Warning: Removed 25932 rows containing missing values (geom_point).

#ggsave(filename=fig.dir %&% "test.tiff",width=3,height=2,units = "in") I can't figure out best parameters
png(filename=fig.dir %&% "Fig-GTEx_TW_EN_CV.png",width=600,height=600)
p
## Warning: Removed 24133 rows containing missing values (geom_point).
## Warning: Removed 24510 rows containing missing values (geom_point).
## Warning: Removed 27268 rows containing missing values (geom_point).
## Warning: Removed 25140 rows containing missing values (geom_point).
## Warning: Removed 25335 rows containing missing values (geom_point).
## Warning: Removed 22325 rows containing missing values (geom_point).
## Warning: Removed 23495 rows containing missing values (geom_point).
## Warning: Removed 23178 rows containing missing values (geom_point).
## Warning: Removed 25932 rows containing missing values (geom_point).
dev.off()
## quartz_off_screen 
##                 2
tiff(filename=fig.dir %&% "Fig-GTEx_TW_EN_CV.tiff",width=600,height=600)
p
## Warning: Removed 24133 rows containing missing values (geom_point).
## Warning: Removed 24510 rows containing missing values (geom_point).
## Warning: Removed 27268 rows containing missing values (geom_point).
## Warning: Removed 25140 rows containing missing values (geom_point).
## Warning: Removed 25335 rows containing missing values (geom_point).
## Warning: Removed 22325 rows containing missing values (geom_point).
## Warning: Removed 23495 rows containing missing values (geom_point).
## Warning: Removed 23178 rows containing missing values (geom_point).
## Warning: Removed 25932 rows containing missing values (geom_point).
dev.off()
## quartz_off_screen 
##                 2

FigS3

GTEx cross-tissue vs tissue-specific h2 (A) and SE (B).

h2TS<-read.table(my.dir %&% "GTEx_Tissue-Specific_local_h2.txt",header=T)
seTS<-read.table(my.dir %&% "GTEx_Tissue-Specific_local_se.txt",header=T)

gh2_TS<-gather(h2TS,"CrossTissue","Tissue",3:11)
colnames(gh2_TS) <- c('ensid','CrossTissue','TissueName','Tissue')
figS2a<-ggplot(gh2_TS,aes(x=Tissue,y=CrossTissue)) +facet_wrap(~TissueName,scales="fixed",ncol=3) + geom_point(alpha=1/10) + geom_abline(intercept=0, slope=1,color='red')  + ylab(expression('Cross-Tissue h'^2)) + xlab(expression('Tissue-Specific h'^2)) +  ggtitle("A")+ coord_cartesian(ylim=c(-0.1,1.1),xlim=c(-0.1,1.1)) + theme(plot.title = element_text(hjust = 0,face="bold"))+theme_bw()

gse_TS<-gather(seTS,"CrossTissue","Tissue",3:11)
colnames(gse_TS) <- c('ensid','CrossTissue','TissueName','Tissue')
figS2b<-ggplot(gse_TS,aes(x=Tissue,y=CrossTissue)) +facet_wrap(~TissueName,scales="fixed",ncol=3) + geom_point(alpha=1/10) + geom_abline(intercept=0, slope=1,color='red') + ylab('Cross-Tissue SE') + xlab('Tissue-Specific SE') +  ggtitle("B") + coord_cartesian(ylim=c(-0.01,0.16),xlim=c(-0.01,0.16))+ theme(plot.title = element_text(hjust = 0,face="bold"))+theme_bw() 
multiplot(figS2a,figS2b)

multiplot(figS2a,figS2b)

tiff(filename=fig.dir %&% "Fig-GTEx-CT-v-TS.tiff",width=480,height=960)
multiplot(figS2a,figS2b)
dev.off()
## quartz_off_screen 
##                 2
png(filename=fig.dir %&% "Fig-GTEx-CT-v-TS.png",width=480,height=960)
multiplot(figS2a,figS2b)
dev.off()
## quartz_off_screen 
##                 2

FigS4

check tissue-specific correlations PGE vs PVE

hun <- read.table(my.vol %&% 'cross-tissue_exp_BSLMM-s100K_iterations_all_chr1-22_2015-07-20.txt',header=T)
ct <- hun %>% arrange(pge50) %>% mutate(position=1:length(pge50),`pge025>0.01`=pge025>0.01) 
ct <- ct[complete.cases(ct),]
h2.TW <- read.table(my.vol %&% 'GTEx_Tissue-Wide_local_h2_se_geneinfo.txt',header=TRUE, sep='\t')
genenames <- h2.TW %>% select(gene=EnsemblGeneID,genename=AssociatedGeneName)
ct <- left_join(ct,genenames,by='gene')

tislist <- scan('/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/nine.tissue.list',sep="\n",what="character")
ts <- data.frame()
for(tis in tislist){
  data <- read.table(my.vol %&% tis %&% '_TS_exp_BSLMM-s100K_iterations_all_chr1-22_2015-08-06.txt',header=T,sep="\t")  
  subdata <- select(data,gene,pve50,pge50,pge025,pge975) %>% mutate(tissue=tis,`pge025>0.01`=pge025>0.01)
  res<-cor.test(subdata$pge50,subdata$pve50)
  #cat(tis,"\tPearson R=",round(res$estimate,3),"\tP-value=",res$p.value,"\n")
  ts <- rbind(ts,subdata)
}

##combine CT and TS PGE vs. PVE in one plot
subct <- select(ct,gene,pve50,pge50,pge025,pge975) %>% mutate(tissue="cross-tissue",`pge025>0.01`=pge025>0.01)
ctts <- rbind(subct,ts) %>% mutate(tissue=factor(tissue,levels=c("cross-tissue","Adipose-Subcutaneous","Artery-Tibial","Heart-LeftVentricle","Lung","Muscle-Skeletal","Nerve-Tibial","Skin-SunExposed(Lowerleg)","Thyroid","WholeBlood")))
p<-ggplot(ctts,aes(x=pve50,y=pge50,ymin=pge025,ymax=pge975,color=`pge025>0.01`)) + facet_wrap(~tissue,ncol=2) + geom_pointrange(col='gray')+geom_point()+theme_bw()+ coord_cartesian(xlim=c(-0.05,1.05))

###calc % nonzero for each tissue
### ADD to plot legend
pvec<-vector()
h2vec<-vector()
a<-ctts %>% select(tissue,`pge025>0.01`) %>% spread(tissue,`pge025>0.01`)
for(i in 1:10){
  tis<-colnames(a)[i]
  per <- signif(table(a[,i])/sum(table(a[,i])),3)*100
  pvec <- c(pvec,per[2])
}
###calc mean PVE for each tissue
a<-ctts %>% select(tissue,pve50) %>% spread(tissue,pve50)
for(i in 1:10){
  tis<-colnames(a)[i]
  meanh2 <- signif(mean(a[,i],na.rm=TRUE),3)
  h2vec <- c(h2vec,meanh2)
}
pvec<-ifelse(is.na(pvec),0,pvec)
ann_text <- data.frame( pge50 = rep(0.08,10), pve50 = rep(0.55,10), percent= pvec, mean_PVE = h2vec, tissue = factor(colnames(a)), pge025=rep(0.9,10), pge975=rep(0.9,10), nonzeroCI=rep(NA,10), se=rep(0.9,10))
p2<-p+geom_text(data=ann_text,aes(label=paste("mean_PVE ==",mean_PVE,sep="")),color="black",show_guide=F,parse=T,hjust=0,size=3) 
ann_text <- data.frame( pge50 = rep(0.20,10), pve50 = rep(0.55,10), percent= pvec, mean_PVE = h2vec, tissue = factor(colnames(a)), pge025=rep(0.9,10), pge975=rep(0.9,10), nonzeroCI=rep(NA,10), se=rep(0.9,10))
p3<-p2+geom_text(data=ann_text,aes(label=paste("percent_TRUE ==",percent,sep="")),color="black",show_guide=F,parse=T,hjust=0,size=3)+ theme(legend.justification=c(1,1), legend.position=c(1,1))
p3

png(filename=fig.dir %&% "Fig-GTEx_CT-TS_BSLMM.png",width=600,height=600)
p3
dev.off()
## quartz_off_screen 
##                 2
tiff(filename=fig.dir %&% "Fig-GTEx_CT-TS_BSLMM.tiff",width=600,height=600)
p3
dev.off()
## quartz_off_screen 
##                 2

FigS5

GTEx cross-tissue and tissue-wide h2 (A) and SE (B).

h2TW<-read.table(my.dir %&% "GTEx_Tissue-Wide_local_h2.txt",header=T)
seTW<-read.table(my.dir %&% "GTEx_Tissue-Wide_local_se.txt",header=T)

gh2_TW<-gather(h2TW,"CrossTissue","Tissue",3:11)
colnames(gh2_TW) <- c('ensid','CrossTissue','TissueName','Tissue')
fig5a<-ggplot(gh2_TW,aes(x=Tissue,y=CrossTissue)) +facet_wrap(~TissueName,scales="fixed",ncol=3) + geom_point(alpha=1/10) + geom_abline(intercept=0, slope=1,color='red')  + ylab(expression('Cross-Tissue h'^2)) + xlab(expression('Tissue-Wide h'^2)) +  ggtitle("A")+ coord_cartesian(ylim=c(-0.1,1.1),xlim=c(-0.1,1.1)) +theme_bw()+ theme(plot.title = element_text(hjust = 0,face="bold")) 

gse_TW<-gather(seTW,"CrossTissue","Tissue",3:11)
colnames(gse_TW) <- c('ensid','CrossTissue','TissueName','Tissue')
fig5b<-ggplot(gse_TW,aes(x=Tissue,y=CrossTissue)) +facet_wrap(~TissueName,scales="fixed",ncol=3) + geom_point(alpha=1/10) + geom_abline(intercept=0, slope=1,color='red')  + ylab('Cross-Tissue SE') + xlab('Tissue-Wide SE') +  ggtitle("B") + coord_cartesian(ylim=c(-0.01,0.16),xlim=c(-0.01,0.16))+theme_bw()+ theme(plot.title = element_text(hjust = 0,face="bold")) 
multiplot(fig5a,fig5b)

tiff(filename=fig.dir %&% "Fig-GTEx-CT-v-TW.tiff",width=480,height=960)
multiplot(fig5a,fig5b)
dev.off()
## quartz_off_screen 
##                 2
png(filename=fig.dir %&% "Fig-GTEx-CT-v-TW.png",width=480,height=960)
multiplot(fig5a,fig5b)
dev.off()
## quartz_off_screen 
##                 2

FigS6

check tissue-wide correlations PGE v PVE

tislist <- scan('/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/nine.tissue.list',sep="\n",what="character")
tw <- data.frame()
for(tis in tislist){
  data <- read.table(my.vol %&% tis %&% '_TW_exp_BSLMM-s100K_iterations_all_chr1-22_2015-10-18.txt',header=T,sep="\t") 
  explist <- scan(out.dir %&% tis %&% ".meanRPKMgt0.1_3samplesRPKMgt0_genelist","c")
  data <- dplyr::filter(data,gene %in% explist)
  subdata <- select(data,gene,pve50,pge50,pge025,pge975) %>% mutate(tissue=tis,`pge025>0.01`=pge025>0.01)
  res<-cor.test(subdata$pge50,subdata$pve50)
  #cat(tis,"\tPearson R=",round(res$estimate,3),"\tP-value=",res$p.value,"\n")
  tw <- rbind(tw,subdata)
}

p<-ggplot(tw,aes(x=pve50,y=pge50,ymin=pge025,ymax=pge975,color=`pge025>0.01`) ) + facet_wrap(~tissue,ncol=3) + geom_pointrange(col='gray')+geom_point()+theme_bw()+ coord_cartesian(xlim=c(-0.05,1.05))

###calc % nonzero for each tissue
### ADD to plot legend
pvec<-vector()
h2vec<-vector()
a<-tw %>% select(tissue,`pge025>0.01`) %>% spread(tissue,`pge025>0.01`)
for(i in 1:9){
  tis<-colnames(a)[i]
  #cat("\n\n---",tis,"---\n")
  #print(table(a[,i]))
  per <- signif(table(a[,i])/sum(table(a[,i])),3)*100
  #print(per)
  pvec <- c(pvec,per[2])
}
###calc mean PVE for each tissue
a<-tw %>% select(tissue,pve50) %>% spread(tissue,pve50)
for(i in 1:9){
  tis<-colnames(a)[i]
  meanh2 <- signif(mean(a[,i],na.rm=TRUE),3)
  #cat(tis,"mean PVE:",meanh2,"\n")
  h2vec <- c(h2vec,meanh2)
}
pvec<-ifelse(is.na(pvec),0,pvec)
ann_text <- data.frame( pge50 = rep(0.08,9), pve50 = rep(0.05,9), percent= pvec, mean_PVE = h2vec, tissue = factor(colnames(a)), pge025=rep(0.9,9), pge975=rep(0.9,9), nonzeroCI=rep(NA,9), se=rep(0.9,9))
p2<-p+geom_text(data=ann_text,aes(label=paste("mean_PVE ==",mean_PVE,sep="")),color="black",show_guide=F,parse=T,hjust=0,size=3) 
ann_text <- data.frame( pge50 = rep(0.20,9), pve50 = rep(0.05,9), percent= pvec, mean_PVE = h2vec, tissue = factor(colnames(a)), pge025=rep(0.9,9), pge975=rep(0.9,9), nonzeroCI=rep(NA,9), se=rep(0.9,9))
p3<-p2+geom_text(data=ann_text,aes(label=paste("percent_TRUE ==",percent,sep="")),color="black",show_guide=F,parse=T,hjust=0,size=3)+ theme(legend.justification=c(1,0), legend.position=c(1,0))
p3

png(filename=fig.dir %&% "Fig-GTEx_TW_BSLMM.png",width=720,height=480)
p3
dev.off()
## quartz_off_screen 
##                 2
tiff(filename=fig.dir %&% "Fig-GTEx_TW_BSLMM.tiff",width=720,height=480)
p3
dev.off()
## quartz_off_screen 
##                 2

FigS7

Plot CT and TS together

cten<-read.table(otd.dir %&% 'cross-tissue_exp_10-foldCV_elasticNet_R2_for_ggplot2.txt',header=T,check.names=F)
ngenesall <- length(unique(cten$gene))
g_cten<-gather(cten,alpha,R2,3:5)

tislist <- scan('/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/nine.tissue.list','c')
finalgdata <- data.frame()
for(tis in tislist){
  alpha1 <- read.table(otd.dir %&% 'TS_' %&% tis %&% '_exp_10-foldCV_elasticNet_alpha1_hapmapSnpsCEU_all_chr1-22_2015-08-27.txt',header=TRUE) %>%  mutate(`1`=R2) %>% select(gene,`1`)
  ngenesall <- length(unique(alpha1$gene))
  alpha95 <- read.table(otd.dir %&% 'TS_' %&% tis %&% '_exp_10-foldCV_elasticNet_alpha0.95_hapmapSnpsCEU_all_chr1-22_2015-08-27.txt',header=TRUE) %>% mutate(`0.95`=R2) %>% select(gene,`0.95`)
  alpha50 <- read.table(otd.dir %&% 'TS_' %&% tis %&% '_exp_10-foldCV_elasticNet_alpha0.5_hapmapSnpsCEU_all_chr1-22_2015-08-27.txt',header=TRUE) %>% mutate(`0.50`=R2) %>% select(gene,`0.50`)
  alpha05 <- read.table(otd.dir %&% 'TS_' %&% tis %&% '_exp_10-foldCV_elasticNet_alpha0.05_hapmapSnpsCEU_all_chr1-22_2015-08-27.txt',header=TRUE) %>% mutate(`0.05`=R2) %>% select(gene,`0.05`)
  data <- inner_join(alpha05,alpha50,by='gene')
  data <- inner_join(data,alpha95,by='gene')
  data <- inner_join(data,alpha1,by='gene')
  gdata <- gather(data,alpha,R2,2:4) %>% mutate(tissue=tis)
  finalgdata <- rbind(finalgdata,gdata)
}

ctsort <- select(g_cten,gene,`1`,alpha,R2) %>% mutate(tissue=g_cten$`cross-tissue`)
ctts <- rbind(ctsort,finalgdata)
p<-ggplot(ctts, aes(y = `1` - R2, x = `1`, group=alpha, color=alpha)) + facet_wrap(~tissue,nrow=2) + geom_point(show_guide = TRUE) + ylab(expression(paste("R"^2, " difference (LASSO - alpha)"))) + xlab(expression(paste("R"^2, " (LASSO)"))) +theme_bw(15)+ theme(legend.justification=c(1,1), legend.position=c(1,1))
p
## Warning: Removed 10623 rows containing missing values (geom_point).
## Warning: Removed 26123 rows containing missing values (geom_point).
## Warning: Removed 24646 rows containing missing values (geom_point).
## Warning: Removed 26515 rows containing missing values (geom_point).
## Warning: Removed 26421 rows containing missing values (geom_point).
## Warning: Removed 25074 rows containing missing values (geom_point).
## Warning: Removed 24244 rows containing missing values (geom_point).
## Warning: Removed 22490 rows containing missing values (geom_point).
## Warning: Removed 23572 rows containing missing values (geom_point).
## Warning: Removed 22612 rows containing missing values (geom_point).

png(filename=fig.dir %&% "Fig-GTEx_CT-TS_EN_CV.png",width=960,height=480)
p
## Warning: Removed 10623 rows containing missing values (geom_point).
## Warning: Removed 26123 rows containing missing values (geom_point).
## Warning: Removed 24646 rows containing missing values (geom_point).
## Warning: Removed 26515 rows containing missing values (geom_point).
## Warning: Removed 26421 rows containing missing values (geom_point).
## Warning: Removed 25074 rows containing missing values (geom_point).
## Warning: Removed 24244 rows containing missing values (geom_point).
## Warning: Removed 22490 rows containing missing values (geom_point).
## Warning: Removed 23572 rows containing missing values (geom_point).
## Warning: Removed 22612 rows containing missing values (geom_point).
dev.off()
## quartz_off_screen 
##                 2
tiff(filename=fig.dir %&% "Fig-GTEx_CT-TS_EN_CV.tiff",width=960,height=480)
p
## Warning: Removed 10623 rows containing missing values (geom_point).
## Warning: Removed 26123 rows containing missing values (geom_point).
## Warning: Removed 24646 rows containing missing values (geom_point).
## Warning: Removed 26515 rows containing missing values (geom_point).
## Warning: Removed 26421 rows containing missing values (geom_point).
## Warning: Removed 25074 rows containing missing values (geom_point).
## Warning: Removed 24244 rows containing missing values (geom_point).
## Warning: Removed 22490 rows containing missing values (geom_point).
## Warning: Removed 23572 rows containing missing values (geom_point).
## Warning: Removed 22612 rows containing missing values (geom_point).
dev.off()
## quartz_off_screen 
##                 2